Exercise

In all of the foregoing, we have assumed a fixed value of the dispersion parameter \(k\), of the negative binomial measurement model. Now relax this assumption: estimate \(k\) along with the other parameters. How much is the fit improved? How do the MLE values of the other parameters change?


Solution

Setup

First, we set up the problem, loading packages, initializing a parallel environment, constructing our SIR model, and identifying the putative MLE. The pomp constructed by the source command below is named measSIR. We take for our putative MLE the point in our database with the highest measured likelihood.

library(pomp)
library(tidyverse)
library(doParallel)
library(doRNG)
registerDoParallel()

source("https://kingaa.github.io/sbied/pfilter/model.R")

read_csv("measles_params.csv") %>%
  filter(!is.na(loglik), loglik.se<1) %>%
  filter(loglik==max(loglik)) %>%
  select(-loglik,-loglik.se) -> coef(measSIR)

The putative MLE is

Beta 3.79e+00
mu_IR 2.00e+00
rho 5.79e-02
k 1.00e+01
eta 5.88e-01
N 3.80e+04

A profile calculation

Because this ridge appears flat in the \(\eta\)-direction, we will profile over \(\eta\) to improve our estimates. The following sets up a design of guesses, to be used as starting points for the profile computation.

read_csv("fitall_params.csv") %>%
  filter(loglik>max(loglik)-10,loglik.se<1) %>%
  sapply(range) -> box

freeze(
  profile_design(
    eta=seq(box[1,"eta"],box[2,"eta"],length=20),
    lower=box[1,c("Beta","rho","k","mu_IR")],
    upper=box[2,c("Beta","rho","k","mu_IR")],
    nprof=15, type="runif"
  ),
  seed=1893696051
)-> guesses
registerDoRNG(830007657)
foreach(guess=iter(guesses,"row"), .combine=rbind) %dopar% {
  library(pomp)
  library(tidyverse)
  mf1 %>%
    mif2(
      Nmif=100,Np=2000,
      params=c(unlist(guess),fixed_params),
      rw.sd=rw.sd(Beta=0.02, rho=0.02, mu_IR=0.02, k=0.02),
      partrans=parameter_trans(log=c("Beta","mu_IR","k"),logit="rho"),
      paramnames=c("Beta","mu_IR","k","rho")
    ) %>%
    mif2(Nmif=100,cooling.fraction.50=0.3) -> mf
  replicate(
    10,
    mf %>% pfilter(Np=2000) %>% logLik()
  ) %>% logmeanexp(se=TRUE) -> ll
  mf %>% coef() %>% bind_rows() %>%
    bind_cols(loglik=ll[1],loglik.se=ll[2])
} -> results
results %>%
  filter(
    is.finite(loglik),
    loglik.se<1
  ) %>%
  group_by(eta) %>%
  filter(rank(-loglik)<=2) %>%
  ungroup() %>%
  gather(variable,value,loglik,mu_IR,rho,Beta,k) %>%
  ggplot(aes(x=eta,y=value))+
  geom_point()+
  labs(y=NULL)+
  facet_wrap(~variable,scales="free_y")


Licensed under the Creative Commons Attribution-NonCommercial license. Please share and remix noncommercially, mentioning its origin.
CC-BY_NC

R codes
Back to the lesson
Course homepage